rho = 0.8
y = c(0, 0)
Niter = 1000
theta = array(0, c(2, Niter))
theta[, 1] = c(5, 8) # starting point, change it to different values
for(k in 2:Niter){
theta[1, k] = y[1] + rho * (theta[2, k-1] - y[2]) + sqrt(1 - rho^2) * rnorm(1)
theta[2, k] = y[2] + rho * (theta[1, k] - y[1]) + sqrt(1 - rho^2) * rnorm(1)
}
# visualize trace (2dim)
plot(t(theta), type = 'l', xlab = expression(theta[1]), ylab = expression(theta[2]))
points(theta[1, 1], theta[2, 1], col = "red", cex = 3)
library(coda)
# visualize each dim
par(mfrow = c(2, 2))
for(k in 1:2){
plot(theta[k, ], type = 'l', xlab = "iteration", ylab = expression(theta), main = paste("gibbs samples for theta", k))
}
for(k in 1:2){
plot(theta[k, 1:30], type = 'l', xlab = "iteration", ylab = expression(theta), main = paste("first 30 iterations, theta", k))
}
# choose burnin
par(mfrow = c(1, 3))
burnin <- 30
theta <- theta[, burnin:dim(theta)[2]]
for(k in 1:2){
plot(theta[k, ], type = 'l', xlab = "iteration", ylab = expression(theta), main = paste("gibbs samples for theta", k, "(after burnin)"))
hist(theta[k, ], xlab = expression(theta), freq = FALSE, main = paste("theta", k))
lines(seq(-3, 3, length.out = 100), dnorm(seq(-3, 3, length.out = 100)), col = "red")
acf(theta[k, ])
}
for(k in 1:2){
print(paste("Effective sample size for theta", k, "is", effectiveSize(theta[k, ])))
}
## [1] "Effective sample size for theta 1 is 209.924170799959"
## [1] "Effective sample size for theta 2 is 220.527568827683"
n = 20
alpha = 0.5
beta = 0.5
Niter = 5000
x = rep(1, Niter)
y = rep(0.5, Niter)
for(k in 2:Niter){
x[k] = rbinom(1, n, y[k-1])
y[k] = rbeta(1, alpha + x[k], n - x[k] + beta)
}
plot(jitter(x), y, pch = 19, cex = 0.2, xlab = 'x', ylab = 'y')
# visualize each dim
par(mfrow = c(2, 2))
plot(y, type = 'l', xlab = "iteration", ylab = "y", main = "gibbs samples for y")
plot(y[1:30], type = 'l', xlab = "iteration", ylab = "y", main = "first 30 iterations")
plot(y[1:80], type = 'l', xlab = "iteration", ylab = "y", main = "first 80 iterations")
plot(y[1:200], type = 'l', xlab = "iteration", ylab = "y", main = "first 200 iterations")
par(mfrow = c(1, 3))
acf(y, main = "no thinning")
acf(y[seq(1, length(y), by = 5)], main = "thin by 5")
acf(y[seq(1, length(y), by = 20)], main = 'thin by 20')
print(paste("Effective sample size for y is", effectiveSize(y)))
## [1] "Effective sample size for y is 139.518663451208"
library(mvtnorm)
mcmc_MH <- function(Niter, proposal_sigma, mu_init, LL){
mu_new <- mu_init
d <- length(mu_new)
samples_MH <- mu_new
accept_MH <- 0
for(iter in 1: Niter){
if(d > 1){
propose_MH_mu <- as.numeric(rmvnorm(1, mean = mu_new, sigma = diag(rep(proposal_sigma^2, d))))
}
if(d == 1){
propose_MH_mu <- rnorm(1, mean = mu_new, sd = proposal_sigma)
}
rtemp <- LL(propose_MH_mu) - LL(mu_new)
if(runif(1) < exp(rtemp)){
mu_new <- propose_MH_mu
accept_MH <- accept_MH + 1
}
samples_MH <- rbind(samples_MH, mu_new)
}
return(list(samples = samples_MH, acceptance = accept_MH /Niter))
}
LL <- function(x){
sum(dnorm(x, log = TRUE))
}
Niter = 1000
proposal_sigma = 2 # try other values
mu_init = c(3, 3) # try other values
samplesMH = mcmc_MH (Niter, proposal_sigma, mu_init, LL)
plot(apply(samplesMH$samples, 2, jitter), type = 'l', xlab = '', ylab = '')
points(t(mu_init), cex = 2, col = 'red')
par(mfrow = c(2, 3))
for(k in 1:2){
plot(jitter(samplesMH$samples[, k]), type = 'l', xlab = '', ylab = '')
hist(samplesMH$samples[, k], xlab = '', main = '', freq = FALSE)
lines(seq(-3, 3, length.out = 100), dnorm(seq(-3, 3, length.out = 100)), col = "red")
acf(samplesMH$samples[, k])
print(paste("Effective sample size", k, " is", effectiveSize(samplesMH$samples[, k])))
}
## [1] "Effective sample size 1 is 126.447072390848"
## [1] "Effective sample size 2 is 170.754497054289"
print(paste("Acceptance rate is", samplesMH$acceptance))
## [1] "Acceptance rate is 0.32"
Niter = 1000
proposal_sigma = 2
mu_init_vals = cbind(c(6, 6, -6, -6), c(6, -6, 6, -6))
burnin = 100
samplesMHmulti = apply(mu_init_vals, 1, function(mu_init) mcmc_MH (Niter, proposal_sigma, mu_init, LL))
idxburn <- burnin:dim(samplesMHmulti[[1]]$samples)[1]
for(jj in 1:4){
par(mfrow = c(1, 1))
samplesMH = samplesMHmulti[[jj]]
plot(apply(samplesMH$samples, 2, jitter), type = 'l', xlab = '', ylab = '', main = paste("Acceptance rate is", samplesMH$acceptance))
points(t(samplesMH$samples[, 1]), cex = 2, col = 'red')
par(mfrow = c(2, 2))
for(k in 1:2){
plot(jitter(samplesMH$samples[, k]), type = 'l', xlab = '', ylab = '', main = paste("dim", k))
plot(jitter(samplesMH$samples[idxburn, k]), type = 'l', xlab = '', ylab = '', main = paste("Effective sample size", k, " is", round(effectiveSize(samplesMH$samples[idxburn, k]), 2)))
hist(samplesMH$samples[idxburn, k], xlab = '', main = '', freq = FALSE)
lines(seq(-3, 3, length.out = 100), dnorm(seq(-3, 3, length.out = 100)), col = "red")
acf(samplesMH$samples[idxburn, k])
}
}
## Geweke diagnostics
for(k in 1:4){
print(geweke.diag(samplesMHmulti[[k]]$samples))
}
##
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5
##
## var1 var2
## 1.496 1.101
##
##
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5
##
## var1 var2
## 1.3886 -0.6951
##
##
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5
##
## var1 var2
## -1.2186 0.4988
##
##
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5
##
## var1 var2
## -0.6942 -1.2734
# after burnin
idxburn <- burnin:dim(samplesMHmulti[[1]]$samples)[1]
for(k in 1:4){
print(geweke.diag(samplesMHmulti[[k]]$samples[idxburn,]))
}
##
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5
##
## var1 var2
## -0.3984 -0.1541
##
##
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5
##
## var1 var2
## -0.03934 0.02760
##
##
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5
##
## var1 var2
## 0.4089 0.3466
##
##
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5
##
## var1 var2
## 0.4322 0.6646
par(mfrow = c(1, 1))
for(dimselect in 1:2){
samples4chains = c()
for(jj in 1:4){
samples4chains <- cbind(samples4chains, samplesMHmulti[[jj]]$samples[, dimselect])
}
plot(samples4chains[, 1], xlab = 'iteration', ylab = '', main = 'compare 4 chains', col = 1, type = "l", ylim = range(samples4chains))
for(jj in 2:4){
lines(samples4chains[, jj], col = jj)
}
legend("topright", col = 1:4, paste("chain", 1:4), lty = rep(1, 4))
idxinit <- 1:30
plot(samples4chains[idxinit, 1], xlab = 'iteration', ylab = '', main = 'compare 4 chains', col = 1, type = "l", ylim = range(samples4chains))
for(jj in 2:4){
lines(samples4chains[idxinit, jj], col = jj)
}
legend("topright", col = 1:4, paste("chain", 1:4), lty = rep(1, 4))
}
# analyze variance
# 1. do not break one chain into two
burnin <- 100
for(dimselect in 1:2){
samples4chains = c()
for(jj in 1:4){
idxburn <- burnin:dim(samplesMHmulti[[jj]]$samples)[1]
samples4chains <- cbind(samples4chains, samplesMHmulti[[jj]]$samples[idxburn, dimselect])
}
within_chain_var <- apply(samples4chains, 2, var)
means_chains <- apply(samples4chains, 2, mean)
between_chain_var <- var(means_chains)
print(paste("dimension", dimselect))
print(paste("Between chain variance is", round(between_chain_var, 3)))
print("Within chain variance is")
print(round(within_chain_var, 3))
}
## [1] "dimension 1"
## [1] "Between chain variance is 0.021"
## [1] "Within chain variance is"
## [1] 1.130 0.808 0.872 1.025
## [1] "dimension 2"
## [1] "Between chain variance is 0.001"
## [1] "Within chain variance is"
## [1] 1.178 0.994 1.061 0.969
# 2. break each chain into two
burnin <- 100
for(dimselect in 1:2){
samples4chains = c()
for(jj in 1:4){
idxburn <- burnin:dim(samplesMHmulti[[jj]]$samples)[1]
samples4chains <- cbind(samples4chains, samplesMHmulti[[jj]]$samples[idxburn, dimselect])
}
halfnum <- floor(dim(samples4chains)[1]/2)
samples4chains <- cbind(samples4chains[1:halfnum, ], samples4chains[(halfnum+1):(2*halfnum), ])
within_chain_var <- apply(samples4chains, 2, var)
means_chains <- apply(samples4chains, 2, mean)
between_chain_var <- var(means_chains)
print(paste("dimension", dimselect))
print(paste("Between chain variance is", round(between_chain_var, 3)))
print("Within chain variance is")
print(round(within_chain_var, 3))
}
## [1] "dimension 1"
## [1] "Between chain variance is 0.025"
## [1] "Within chain variance is"
## [1] 1.048 0.755 0.770 1.101 1.211 0.839 0.965 0.945
## [1] "dimension 2"
## [1] "Between chain variance is 0.007"
## [1] "Within chain variance is"
## [1] 1.078 1.010 1.123 1.088 1.240 0.979 1.000 0.851
mcmclist4chains <- mcmc.list(as.mcmc(samplesMHmulti[[1]]$samples), as.mcmc(samplesMHmulti[[2]]$samples), as.mcmc(samplesMHmulti[[3]]$samples), as.mcmc(samplesMHmulti[[4]]$samples))
gelman.diag(mcmclist4chains)
## Potential scale reduction factors:
##
## Point est. Upper C.I.
## [1,] 1.01 1.03
## [2,] 1.00 1.01
##
## Multivariate psrf
##
## 1.01
gelman.plot(mcmclist4chains)
## Geweke diagnostics
for(k in 1:4){
print(geweke.diag(mcmclist4chains[[k]]))
geweke.plot(mcmclist4chains[[k]])
}
##
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5
##
## var1 var2
## 1.496 1.101
##
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5
##
## var1 var2
## 1.3886 -0.6951
##
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5
##
## var1 var2
## -1.2186 0.4988
##
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5
##
## var1 var2
## -0.6942 -1.2734
# after burnin
idxburn <- burnin:dim(samplesMHmulti[[1]]$samples)[1]
for(k in 1:4){
print(geweke.diag(samplesMHmulti[[k]]$samples[idxburn,]))
}
##
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5
##
## var1 var2
## -0.3984 -0.1541
##
##
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5
##
## var1 var2
## -0.03934 0.02760
##
##
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5
##
## var1 var2
## 0.4089 0.3466
##
##
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5
##
## var1 var2
## 0.4322 0.6646
Niter = 1000
proposal_sigma = 2
mu_init = c(0, 0)
burnin = 100
proposal_sigma_vals <- c(0.3, 1, 3, 6)
samplesMHmulti = sapply(proposal_sigma_vals, function(proposal_sigma) mcmc_MH (Niter, proposal_sigma, mu_init, LL), simplify = FALSE)
idxburn <- burnin:dim(samplesMHmulti[[1]]$samples)[1]
for(jj in 1:4){
par(mfrow = c(1, 1))
samplesMH = samplesMHmulti[[jj]]
plot(apply(samplesMH$samples, 2, jitter), type = 'l', xlab = '', ylab = '', main = paste("Acceptance rate is", samplesMH$acceptance))
points(t(samplesMH$samples[, 1]), cex = 2, col = 'red')
par(mfrow = c(2, 2))
for(k in 1:2){
plot(jitter(samplesMH$samples[, k]), type = 'l', xlab = '', ylab = '')
plot(jitter(samplesMH$samples[idxburn, k]), type = 'l', xlab = '', ylab = '', main = paste("Effective sample size", k, " is", round(effectiveSize(samplesMH$samples[idxburn, k]), 2)))
hist(samplesMH$samples[idxburn, k], xlab = '', main = '', freq = FALSE)
lines(seq(-3, 3, length.out = 100), dnorm(seq(-3, 3, length.out = 100)), col = "red")
acf(samplesMH$samples[idxburn, k], main = '')
}
}
par(mfrow = c(1, 1))
for(dimselect in 1:2){
samples4chains = c()
for(jj in 1:4){
samples4chains <- cbind(samples4chains, samplesMHmulti[[jj]]$samples[, dimselect])
}
plot(samples4chains[, 1], xlab = 'iteration', ylab = '', main = 'compare 4 chains', col = 1, type = "l", ylim = range(samples4chains))
for(jj in 2:4){
lines(samples4chains[, jj], col = jj)
}
legend("topright", col = 1:4, paste("chain", 1:4), lty = rep(1, 4))
idxinit <- 1:30
plot(samples4chains[idxinit, 1], xlab = 'iteration', ylab = '', main = 'compare 4 chains', col = 1, type = "l", ylim = range(samples4chains))
for(jj in 2:4){
lines(samples4chains[idxinit, jj], col = jj)
}
legend("topright", col = 1:4, paste("chain", 1:4), lty = rep(1, 4))
}
# Calculate R hat
burnin <- 100
for(dimselect in 1:2){
samples4chains = c()
for(jj in 1:4){
idxburn <- burnin:dim(samplesMHmulti[[jj]]$samples)[1]
samples4chains <- cbind(samples4chains, samplesMHmulti[[jj]]$samples[idxburn, dimselect])
}
halfnum <- floor(dim(samples4chains)[1]/2)
samples4chains <- cbind(samples4chains[1:halfnum, ], samples4chains[(halfnum+1):(2*halfnum), ])
within_chain_var <- apply(samples4chains, 2, var)
means_chains <- apply(samples4chains, 2, mean)
between_chain_var <- var(means_chains)
print(paste("dimension", dimselect))
print(paste("Between chain variance is", round(between_chain_var, 3)))
print("Within chain variance is")
print(round(within_chain_var, 3))
}
## [1] "dimension 1"
## [1] "Between chain variance is 0.134"
## [1] "Within chain variance is"
## [1] 0.559 1.081 0.587 0.322 1.179 0.993 0.669 0.824
## [1] "dimension 2"
## [1] "Between chain variance is 0.011"
## [1] "Within chain variance is"
## [1] 0.703 0.889 0.722 1.057 0.840 0.793 0.939 0.990
mcmclist4chains <- mcmc.list(as.mcmc(samplesMHmulti[[1]]$samples), as.mcmc(samplesMHmulti[[2]]$samples), as.mcmc(samplesMHmulti[[3]]$samples), as.mcmc(samplesMHmulti[[4]]$samples))
gelman.diag(mcmclist4chains)
## Potential scale reduction factors:
##
## Point est. Upper C.I.
## [1,] 1.11 1.29
## [2,] 1.01 1.05
##
## Multivariate psrf
##
## 1.12
gelman.plot(mcmclist4chains)
## Geweke diagnostics
for(k in 1:4){
print(geweke.diag(mcmclist4chains[[k]]))
geweke.plot(mcmclist4chains[[k]])
}
##
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5
##
## var1 var2
## -0.8755 -2.0223
##
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5
##
## var1 var2
## 1.1565 -0.6244
##
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5
##
## var1 var2
## -0.6528 0.2308
##
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5
##
## var1 var2
## -0.09499 -0.27813
# after burnin
idxburn <- burnin:dim(samplesMHmulti[[1]]$samples)[1]
for(k in 1:4){
print(geweke.diag(samplesMHmulti[[k]]$samples[idxburn,]))
}
##
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5
##
## var1 var2
## -3.4830 0.6443
##
##
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5
##
## var1 var2
## -0.6974 -0.1809
##
##
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5
##
## var1 var2
## 0.3939 1.4149
##
##
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5
##
## var1 var2
## 0.8648 2.1274
Niter = 1000
proposal_sigma = 2
mu_init = c(0, 0)
burnin = 100
samplesMHmulti = sapply(1:4, function(proposal_sigma) mcmc_MH (Niter, proposal_sigma, mu_init, LL), simplify = FALSE)
# 1. do not break one chain into two
burnin <- 100
idxburn <- burnin:dim(samplesMHmulti[[1]]$samples)[1]
for(dimselect in 1:2){
samples4chains = c()
for(jj in 1:4){
idxburn <- burnin:dim(samplesMHmulti[[jj]]$samples)[1]
samples4chains <- cbind(samples4chains, samplesMHmulti[[jj]]$samples[idxburn, dimselect])
}
within_chain_var <- mean(apply(samples4chains, 2, var))
means_chains <- apply(samples4chains, 2, mean)
n = length(idxburn)
between_chain_var <- var(means_chains) * n
Rstat <- sqrt(((1 - 1/n) * within_chain_var + between_chain_var /n)/within_chain_var)
print(paste("dimension", dimselect))
print(paste("Between chain variance is", round(between_chain_var, 3)))
print(paste("Within chain variance is", round(within_chain_var, 3)))
print(paste("Rhat is", round(Rstat, 3)))
}
## [1] "dimension 1"
## [1] "Between chain variance is 7.804"
## [1] "Within chain variance is 0.886"
## [1] "Rhat is 1.004"
## [1] "dimension 2"
## [1] "Between chain variance is 13.743"
## [1] "Within chain variance is 0.943"
## [1] "Rhat is 1.007"
# 2. break each chain into two
burnin <- 100
for(dimselect in 1:2){
samples4chains = c()
for(jj in 1:4){
idxburn <- burnin:dim(samplesMHmulti[[jj]]$samples)[1]
samples4chains <- cbind(samples4chains, samplesMHmulti[[jj]]$samples[idxburn, dimselect])
}
halfnum <- floor(dim(samples4chains)[1]/2)
samples4chains <- cbind(samples4chains[1:halfnum, ], samples4chains[(halfnum+1):(2*halfnum), ])
within_chain_var <- mean(apply(samples4chains, 2, var))
means_chains <- apply(samples4chains, 2, mean)
n = halfnum
between_chain_var <- var(means_chains) * n
Rstat <- sqrt(((1 - 1/n) * within_chain_var + between_chain_var / n)/within_chain_var)
print(paste("dimension", dimselect))
print(paste("Between chain variance is", round(between_chain_var, 3)))
print(paste("Within chain variance is", round(within_chain_var, 3)))
print(paste("Rhat is", round(Rstat, 3)))
}
## [1] "dimension 1"
## [1] "Between chain variance is 6.411"
## [1] "Within chain variance is 0.881"
## [1] "Rhat is 1.007"
## [1] "dimension 2"
## [1] "Between chain variance is 10.625"
## [1] "Within chain variance is 0.935"
## [1] "Rhat is 1.011"
mcmclist4chains <- mcmc.list(as.mcmc(samplesMHmulti[[1]]$samples), as.mcmc(samplesMHmulti[[2]]$samples), as.mcmc(samplesMHmulti[[3]]$samples), as.mcmc(samplesMHmulti[[4]]$samples))
gelman.diag(mcmclist4chains)
## Potential scale reduction factors:
##
## Point est. Upper C.I.
## [1,] 1.01 1.04
## [2,] 1.02 1.05
##
## Multivariate psrf
##
## 1.02
gelman.plot(mcmclist4chains)
## Geweke diagnostics
for(k in 1:4){
print(geweke.diag(mcmclist4chains[[k]]))
geweke.plot(mcmclist4chains[[k]])
}
##
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5
##
## var1 var2
## 0.07346 -0.82772
##
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5
##
## var1 var2
## 0.7323 -0.2362
##
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5
##
## var1 var2
## -0.9449 0.5986
##
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5
##
## var1 var2
## 0.5923 -0.7085
# after burnin
idxburn <- burnin:dim(samplesMHmulti[[1]]$samples)[1]
for(k in 1:4){
print(geweke.diag(samplesMHmulti[[k]]$samples[idxburn,]))
}
##
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5
##
## var1 var2
## -0.2862 -0.7265
##
##
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5
##
## var1 var2
## 0.1408 -0.9675
##
##
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5
##
## var1 var2
## -0.1198 -0.9092
##
##
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5
##
## var1 var2
## -1.4935 0.1664